FULLY ADAPTIVE MULTIRESOLUTION SCHEMES FOR STRONGLY 
DEGENERATE PARABOLIC EQUATIONS WITH DISCONTINUOUS FLUX 
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Abstract. A fully adaptive finite volume mult iresolut ion scheme for one-dimensional strongly degenerate 

parabolic equations with discontinuous flux is presented. The numerical scheme is based on a finite volume 

^«s discretization using the approximation of Engquist-Osher for the flux and explicit time stepping. An adaptive 

/"— >> multiresolution scheme with cell averages is then used to speed up CPU time and memory requirements. 

»»^ A particular feature of our scheme is the storage of the multiresolution representation of the solution in a 

dynamic graded tree, for sake of data compression and to facilitate navigation. Applications to traffic flow 

,_^ with driver reaction and a clarifier-thickener model illustrate the efficiency of this method. 



oo 



< 



> 

X 



1. Introduction 

1.1. Scope of the paper. High resolution finite volume schemes for the approximation of discontinuous 
solutions to conservation laws are of at least second-order accuracy in regions where the solution is smooth 
and resolve discontinuities sharply and without spurious oscillations. Methods of this type include the 



Cd schemes described in [1, 2, 3, 4, 5, 6]. In standard situations, the solution u{x,t) of a conservation law 

£j (1.1) ut + f{u)^ = 0, {x,t)eQT-^nx[0,T], r!CM 

^^ exhibits strong variations (shocks) in small regions but behaves smoothly on the major portion of the com- 

^ putational domain. The multiresolution technique adaptively concentrates computational effort associated 

(^ with a high resolution scheme on the regions of strong variation. It goes back to Harten [ : ] for hyperbolic 

^D equations and was used by Bihari and Harten [8] and Roussel et al. [U] for parabolic equations. Important 

^r contributions to the analysis of multiresolution methods for conservation laws include [10, 11, iz]. 
^^ In this paper, we present a fully adaptive multiresolution scheme and corresponding numerical experiments 

r ■ for strongly degenerate parabolic equations with discontinuous flux. Specifically, we consider equations of 

^^ the type 

O (1.2) ut + f{'r{x),u)^={-fi{x)A{u),,)^ forxenT:=Mx(0,r], 



where we assume that for each a;, the function /(^{x),-) : M ^ M is piecewise smooth and Lipschitz 

S^ continuous, and that 'y{x) is a vector of scalar parameters that are discontinuous at most at a finite number 

5—1 of points. On the other hand, we assume that the integrated diffusion function A{-) is Lipschitz continuous 

and piecewise smooth with A{v) ^ A{u) for v > u. We admit intervals [a, (3] with A{u) = const, for all 

u G [a,f3], such that (1.2) degenerates into the first-order equation 

(1.3) ut + f{j{x),u)^^0 

wherever u € [a, /?]. If degeneracy occurs on w-intervals of positive length (and not only at isolated points), 
(1.2) is called strongly degenerate. Clearly, solutions of (1.2) are in general discontinuous, and need to be 
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characterized as weak solutions along with an entropy condition. Applications of (1.2) with constant param- 
eters include models of sedimentation-consolidation processes of particulate suspensions [l-i, 14], two-phase 
flow in porous media [I'l], and traffic flow with driver reaction [1(>, IT]. Applications with a discontinuous 
parameter vector 'y{x) include models of traffic flow on highways with discontinuous road surface condi- 
tions [IG, IS], and a model of clarifier-thickener units used in engineering applications for the continuous 
solid-liquid separation of suspensions [I'l, '''>]. In the latter application, the function A{u) models sediment 
compressibility; the special case A = 0, in which we fall back to (1.3), corresponds to a so-called ideal 
suspension of rigid spherical particles forming incompressible sediments. See [''•'] for further applications. 

The novelty of the present paper is that we apply an adaptive multiresolution method to one-dimensional 
initial- value problems for (1.2). This equation is discretized in space by a first-order conservative finite volume 
scheme using the Engquist-Osher approximation, for which convergence results for our class of problems are 
available ["i, I'l, 20, 21, 22]. For time discretization an explicit Euler scheme is used. The multiresolution 
representation of the solution allows to introduce a locally refined mesh by thresholding of the wavelet 
coefficients while controlling the error of the approximation. This allows us to reduce the number of costly 
flux evaluations with respect to the finite volume scheme on a regular fine grid. Hence, a gain in CPU time 
can be obtained. Furthermore, the data are efficiently represented in a dynamic graded tree data structure, 
which also leads to memory compression. 



1.2. Multiresolution schemes. In the following, we briefly outline the underlying ideas of multiresolution 
schemes for conservation laws and parabolic equations. The starting point is a conservative high-resolution 
finite volume discretization on a uniform mesh of (1.1) or (1.3). A multiresolution analysis of the solution 
with subsequent thresholding of the coefficients allows an approximation with less coefficients within a given 
tolerance. This allows us to reduce the number of costly flux evaluations required by the high-resolution 
scheme, which results in a gain of efficiency. For this purpose, either point values or cell averages of the 
numerical solution are defined on a hierarchical sequence of nested dyadic grids. Applying a multiresolution 
analysis to the solution, which can be efficiently done using the fast wavelet transform, we can construct a 
truncated representation by simple thresholding of the obtained coefficients. This procedure yields an efficient 
representation of the solution on a locally refined grid while controlling the error of the approximation. 

The principle of multiresolution data representation consists in considering grid averages of the data 
at different resolutions from the finest to the coarsest grid, and in encoding the differences between two 
grids. Finally, one retains only the grid averages on the coarsest grid and the set of errors (or details) for 
predicting the grid averages of each resolution level in this hierarchy from those of the next coarsest one. In 
regions where the solution is sufficiently smooth, the multiresolution coefficients are small and can hence be 
neglected. Thus, data can be compressed by a thresholding or truncation operation, i.e. by setting to zero 
those components of the representation whose multiresolution coefficients (also called wavelet coefficients or 
details) are in absolute value smaller than a prescribed tolerance. Thresholding allows to control the so- 
called perturbation error thanks to norm equivalences. The representation of the solution in physical space 
corresponds to a locally refined grid. 

The multiresolution analysis of the numerical solution automatically detects discontinuities, since a wavelet 
coefficient takes into account the regularity of a function in each position and on each scale. In [7, 23, 24] 
Harten explored this idea and introduced multiresolution schemes for efficiently solving hyperbolic conserva- 
tions laws. Using the multiresolution representation of the solution, he devised a sensor to decide at which 
positions of a fine mesh the flux should be exactly evaluated, and where otherwise it can be obtained more 
cheaply by interpolation of pre-calculated fluxes on coarser scales. Still in the context of hyperbolic con- 
servation laws and preserving flux evaluations for all fine grid positions, Bihari and Harten [>■] developed a 
second-order adaptive switch for flux evaluations, keeping an essentially non-oscillatory (ENO) scheme where 
multiresolution coefficients were larger than a given tolerance, and otherwise using interpolation. In [25], 
Daubechies wavelets were used as a grid refinement strategy associated with finite difference stencils on an 
irregular grid for solving hyperbolic equations. Centered finite differences are used in [2(3] for approximating 
space derivatives on sparse point approximations (SPR) obtained by interpolating wavelet transforms. An 
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SPR-based multiresolution WENO scheme is presented in [ '.]. For parabolic PDEs a finite volume method 
with dynamical adaptation strategy to advance the grid was developed in [ ]. 

An alternative adaption strategy could be based on local a posteriori error estimates by means of residual 
error computation. Results of a posteriori error estimates have been reported in the literature for elliptic 
problems (see [2H]), parabolic problems (see [29], [30]) and hyperboHc problems (see [-'jl], [32] and the 
references therein), but there is not known results for strongly degenerate parabolic problems. In this sense, 
to compute local error estimator is not easy to be realized in practice, and we prefer to concentrate our effort 
in the strategy based on multiscale technique. The multiresolution strategy proposed herein for strongly 
degenerate parabolic equations with a discontinuous flux produces a gain in computational time and in 
memory. The solution is efficiently represented using a graded tree data structure and the costly fluxes 
are computed on the locally refined grid only. The computational efficiency of the multiresolution method 
is related to the data compression rate, that is, to the amount of significant information preserved after 
thresholding in comparison with the number of grid points of the finest mesh. Thus, efficiency is measured 
in terms of the compression rate and CPU time. 

Finally, although we limit our treatment to one space dimension, the multiresolution scheme can be ex- 
tended to higher dimensional problems in different ways. One possibility is to use higher dimensional wavelet 
transforms constructed by a tensor product approach, and through interpolations of the numerical divergence 
in the sense of cell averages from coarser to finer levels, the method of predicting values hierarchically can 
be extended as done in [ ] . Another possibility is to explore the splitting capability of the divergence by 
directions as in [ ] . Fully three-dimensional computations of fiame instabilities are presented in [35] . 

1.3. Strongly degenerate parabolic equations and conservation laws with discontinuous flux. 

Equation (1.2) combines two independent non-standard ingredients of conservation laws: the strongly de- 
generate diffusion term A{u)xx, and the flux f{'j{x),u) that depends discontinuously on the spatial position x. 
We briefly review some recent results for equations that include either ingredient. 
The basic difficulty associated with degenerate parabolic equations of the type 

(1.4) ut + /(u), = A{u),.,, xencR, t e (0, T] 

is that their solutions need to be defined as weak, in general discontinuous solutions along with an entropy 
condition to select the physically relevant weak solution. In [ ' '] the existence of BV entropy weak solutions 
to an initial-boundary value problem for (1.4) in the sense of Kruzkov [ ] and Vol'pert and Hudjaev [37, 38] is 
shown via the vanishing viscosity method, while their uniqueness is shown by a technique due to Carrillo [.'^i!)]. 
The well-posedness of multi-dimensional Dirichlet initial-boundary value problems for strongly degenerate 
parabolic equations is shown in [^i-, i]. Further recent contributions to the analysis of strongly degenerate 
parabolic equations include [41, 42, 43, 44]. 

Evje and Karlsen [ ' ] show that explicit monotone finite difference schemes [4(1] converge to BV entropy 
solutions for the Cauchy problem for (1.4). These results are extended to several space dimensions in [^V]. 
The convergence of finite volume schemes for initial-boundary value problems is proved in [44, 48]. The 
monotone scheme used for numerical experiments in [19, 20] is the robust Engquist-Osher scheme [49]. 
Thus, (1.4) admits a rigorous convergence analysis for suitable numerical schemes. 

In the context of the clarifier-thickener model, the analysis of (1.2) for the case A = 0, that is, of the 
first-order conservation law with discontinuous flux (1.3), has been topic of a recent series of papers including 
[19, 50, 51], in which a rigorous mathematical (existence and uniqueness) and numerical analysis is provided. 
The main ingredient in these clarifier-thickener models is equation (1.3), where the (with respect to u, 
nonconvex) fiux / and the discontinuous vector- valued coefficient 7 = (71,72) are given functions. When 7 
is smooth, Kruzkov's theory [ ] ensures the existence of a unique and stable entropy weak solution to (1.3). 
Kruzkov's theory does not apply when 7 is discontinuous. In [ ], a variant of Kruzkov's notion of entropy 
weak solution for (1.3) that accounts for the discontinuities in 7 is introduced and existence and uniqueness 
(stability) of such entropy solutions in a certain functional class are proved. The existence of such solutions 
follows from the convergence of various numerical schemes such as front tracking [50], a relaxation scheme 
[51, 52], and upwind difference schemes [19]. 
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Strongly degenerate parabolic equations with discontinuous fluxes are studied in [21, 22, 53]. In [21] 
equations like (1.2) are studied with a concave convective flux u h^ f{'y{x),u) and with {-ji^x) A{u)x)^ 
replaced by A(u)a:x- Existence of an entropy weak solution is established by proving convergence of a 
difference scheme of the type discussed in this paper. Uniqueness and stability issues for entropy weak 
solutions are studied in [22] for a particular class of equations. These analyses are extended to the traffic 
and clariflcr-thickener models studied herein in [ ] and [20], respectively. 

1.4. Time discretization, space discretization, and numerical stability. The numerical scheme for 
the solution of (1.2) is described in [ ]. In this work, the basic scheme is flrst order in time and space. 
We utilize a simple explicit Euler discretization in time. The spatial discretization is done by using the 
Engquist-Osher approximation for the convective part of the flux combined with a second-order conservative 
discretization of the diffusion term. For stability we need to satisfy a CFL condition requiring that in general 
At/{Axy be bounded. In some cases without diffusion (Example 2 of Section 6) we need only that At/ Ax 
be bounded. 

1.5. Outline of this paper. The remainder of this paper is organized as follows. In Section 2, we briefly 
outline two applicative models that lead to an equation of the type (1.2), namely, a model of traffic flow 
with driver reaction and discontinuous road surface conditions (Section 2.1) and a clarifier-thickener model 
(Section 2.2). For detailed derivations of both models, we refer to [16] and [58], respectively. In Section 3, 
we describe the basic numerical finite volume scheme for the discretization of (1.2) on a uniform grid. 

In Section 4, the conservative adaptive multiresolution discretization is introduced. Details on the nu- 
merical method and on its implementation using dynamical data structures can be found in [9]. For the 
particular application to strongly degenerate parabolic equations with a flux that depends on u but not on x, 
we refer to [59]. 

The basic motivation of this approach is to accelerate a given finite volume scheme on a uniform grid 
without loosing accuracy. The principle of the multiresolution analysis is to represent a set of data given 
on a fine grid as values on a coarser grid plus a series of differences at different levels of nested dyadic 
grids. These differences contain the information of the solution when going from a coarse to a finer grid. An 
appealing feature of this data representation is that coefficients are small in regions where the solution is 
smooth. Applying a thresholding of small coefficients a locally refined adaptive grid is defined. The threshold 
is chosen in such a way to guarantee that the discretization error of the reference scheme is balanced with 
the accumulated thresholding error which is introduced in each time step. This yields a memory and CPU 
time reduction while controlling the precision of the computations. The dynamic graded tree is introduced in 
Section 4.1, while the multiresolution transform of a function, which is stored in the graded tree, is outlined 
in Section 4.2. The complete multiresolution algorithm is outlined in Section 4.3. 

An error analysis, which has been adapted from Cohen et al. [11] and is also advanced in [ ] for strongly 
degenerate parabolic equations of the type (1.4), is presented in Section 5. This error analysis motivates the 
choice of two parameters in the thresholding algorithm. In Section 6 we present three numerical examples, 
namely the traffic model (Example 1, Section 6.1), a sub-case of the clarifier-thickener model with A = 
that illustrates the application of the method to (1.3) (Example 2, Section 6.2), and the clarifier-thickener 
model treating a flocculated suspension, now again with a degenerate diffusion term A ^ (Example 3, 
Section 6.3). Numerical results, limitations and extensions of the method are discussed in Section 7. 

2. Applications of strongly degenerate parabolic equations 

2.1. Traffic flow with driver reaction and discontinuous road surface conditions. The classical 
Lighthill-Whitham- Richards (LWR) kinematic wave model [54, 55] for unidirectional traffic flow on a single- 
lane highway starts from the principle of "conservation of cars" ut + {uv)x — for x G M and t > 0, where u 
is the density of cars as a function of distance x and time t and v = v{x,t) is the velocity of the car located 
at position x at time t. The decisive constitutive assumption of the LWR model is that -y is a function of 
u only, V = v{u). In other words, it is assumed that each driver instantaneously adjusts his velocity to the 
local car density. A common choice is v{u) = fmaxV^(u), where Wmax is a maximum velocity a driver assumes 
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on a free highway, and V{u) is a hindrance function taking into account the presence of other cars that urges 
each driver to adjust his speed. Thus, the flux is 

(2.1) f{u) := uv{u) ~ v-iaaxuV [u) for ^ M ^ Wmax, /("") = Otherwise, 

where Uniax is the maximum "bumper-to-bumper" car density. The simplest choice is the hnear interpolation 
V{u) = Vi{u) := 1 — u/umax] but we may also consider the alternative Dick-Greenberg model [50, 57] 

(2.2) V{u) = V2{u) := min{l, Cln(u,„ax/u)}, C > 0. 

The diffusively corrected kinematic wave model (DCKWM) [Ifi, 17] extends the LWR model by a strongly 
degenerating diffusion term. This model incorporates a reaction time r, representing drivers' delay in 
their response to events, and an anticipating distance La, which means that drivers adjust their veloc- 
ity to the density seen an anticipating distance La ahead. In fact, we adopt the equation [17] La = 
max{(ti(M))^/(2a), iinin}, where the first argument is the distance required to decelerate to full stop from 
speed v{u) at deceleration a, and the second imposes a minimal anticipation distance, regardless of how 
small the velocity is. If one assumes that the effects of reaction time and anticipation are only relevant when 
the local car density exceeds a critical value Uc, then the final governing equation (replacing ut + f(u)x — 0) 
of the DCKWM is the strongly degenerate parabolic equation 



(2.3) ut + f{u)x = A{u)xx, xeM, t>Q- A{u):^ a{s)ds, 

Jo 

where (see [16, 58] for details of the derivation) 

(2.4) f{u) = v^a^uV{u), a{u):=\^ vu U r ( ^ ^ T/'r ^^ T " 5 "" 

[-MWmaxl/ [Uj[La[u) + TV^i^^uV [u)) it U > U^. 

(A critical density Uc > automatically arises from the use of (2.2); obviously, V2{u) = for u ^ Uc := 
Mmaxexp(-1/C), SO that (2.4) holds for V{u) = V2{u).) 

We assume that V{u) is chosen such that D'{u) > for Uc < w < Umax- Consequently, the right-hand 
side of (2.3) vanishes on the interval [0, Uc], and possibly at the maximum density u^ax- Thus, the governing 
equation of the DCKWM model (2.3) is strongly degenerate parabolic. 

Following Mochon [±6], Burger and Karlsen [ ' '-] extend the DCKWM traffic model to variable road surface 
conditions by replacing the coefficient Wmax in /(«) = v-n-iaxuV {u) by a discontinuously varying function 
f'max = Vma.-x.{x). However, the degenerate diffusion term models driver psychology and should therefore not 
depend on road surface conditions. Consequently, the new model equation for the traffic model is 

(2.5) Ut + !{-i{x),u)^^ A{u).j,j:, !{-i{x),u) ■.= -i{x)uV{u), 7(a;) := w„iax(a;)- 

For simplicity, we assume that on the major part of the highway, the maximum velocity assumes a constant 
value timax; which is also used as the value of t^max entering the definition of A{u) in (2.4), and that there is 
an interval [a, b] on which the maximum velocity assumes an exceptional value w^ax ¥" ''^max- 



(2.6) v^U^) = ^{x) = { - 



for X G [a, 6], 



w^ax otherwise. 



The initial-value problem for equation (2.5) with Cauchy data u{x,Q) = uo(x) for a; € M is well posed 
[Ki], but we here insist on using a finite domain that can completely be represented by our data structure. 
Therefore we consider a circular road of length L, the initial condition 

(2.7) u(a;,0) = uo(a;), 2:e[0,L], 
and the periodic boundary condition 

(2.8) u(0,i) =u(i,i), ie(0,r]. 

Consequently, the "traffic model" is defined by the periodic initial-boundary value problem (2.5), (2.7), (2.8) 
under the assumptions (2.4) and (2.6), where we assume < a < 5 < L. 
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Figure 1. The clarifier-thickener model. 



2.2. Clarifier-thickener model. The analysis of (1.4) has in part been motivated by a theory of sedi- 
mentation-consolidation processes of flocculated suspensions [13, 20], in which the unknown is the solids 
concentration m as a function of time t and depth x. The particular suspension is characterized by the 
hindered settling function f{u) and the integrated diffusion coefficient A(u), which models the sediment 
compressibility. The function f{u) is assumed to be continuous and piecewise smooth with f{u) > for 
M S (0, Wmax) and f{u) = for w ^ and u ^ Umax- A typical example is 





(2.9) 



/(") 



for U € (0,Umax), 

otherwise, 



> 0, C > 0, 



where Vaa > is the settling velocity of a single particle in unbounded fluid. Moreover, we have that 

f{uK{u) 



(2.10) 



A{u) 



a{s) ds, a{u) := 



10 ^g9u 

Here, A^ > is the solid-fluid density difference, g is the acceleration of gravity, and (Tg{u) is the derivative 
of the material specific effective solid stress function adu). We assume that the solid particles touch each 
other at a critical concentration value (or gel point) ^ Uc ^ Umax, and that 



(2.11) 



a^{u),cr^{u) 



= for u ^ Uc 
> for u > Uc 



This implies that a{u) — for u ^ u^, such that also this application motivates a strongly degenerate 
parabolic equation (1.3). A typical function is 



(2.12) 



(Te(u) = 







for u ^ Uc 



aa[{u/uc)^ — 1] for u > 



(To > 0, ^ > 1. 



The extension of the one-dimensional sedimentation-consolidation equation (1.4) (if f{u) and A{u) have 
the interpretation given herein) to continuous sedimentation processes leads to the so-called clarifier-thickener 
model [20], see Figure 1. We consider a cylindrical vessel of constant cross-sectional area S, which occupies 
the depth interval [a;L,a;R] with xl < and sr, > 0. At depth x = 0, fresh suspension of a given feed 
concentration up G [0,Mmax] is pumped into the unit at a volume rate Qf > 0. Within the unit, the feed 
fiow is divided into an upwards directed and a downwards-directed bulk flow with the signed volume rates 
Ql ^ and Qr ^ 0, where conservation of suspension implies Qf — Qr ~ Qh- Furthermore, we assume 
that the feed suspension is loaded with solids at the given feed concentration uf- Finally, at a; = a;L and 
X — x^, overflow and underflow pipes are provided through which the material leaves the clarifier-thickener 
unit. We assume that the solid and the fiuid phases move at the same velocity through these pipes, so that 
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the solid-fluid relative velocity is zero for x < x\^ and x > xr, which means that the term /(u) — A(u)^ is 
"switched off" outside [xljXr]. See [20] for details. 

We only consider vessels with a constant interior cross-sectional area S and define the velocities ql '■— Qi^lS 
and (7r :— Qk/S. Then the final clarifier-thickener model is given by (1.2) with 

(2.13) f{j{x),u)^j2ix)iu~UF) + liix)f{u), 
where we use the discontinuous parameters 

,r,..^ , N jl for a; e (a;L,XR), JqL for a; < 0, 

(2.14) 7i(a^) = <„ ,, . 72(2;) -< , ^^ 

10 otherwise, [qr tor a; > 0. 

We assume the initial concentration distribution 

(2.15) u(x,0) = ^0(2;), 2: e M; "0(2^) G [0, Umax]- 

Thus, the clarifier-thickener model is specified by (1.2) with the discontinuous fiuxes defined by the continuous 
functions u 1— > f{u),A{u) given by (2.9) and (2.10), the discontinuous parameters (2.13), (2.14), and the 
initial condition (2.15). 

3. Numerical scheme 

The numerical scheme for the solution of (1.2) is essentially described in [:?'i]. We begin the definition of 
the base algorithm discretizing M into cells Ij :— [2^^-1/2, 2:^+1/2)1 where 2;j_|.i/2 — (j + l/2)Aa; with j G Z. 
Let A — At/ Ax, ji — At/{AxY and C/" — uq{xj). For n > we define the approximations according to 

(3.1) [/;+! = c/; - AA_;i(7^.+i/2, u^+,,w^) + mA- (7i,,+i/2A+a(c/;)) , 

where 

(3-2) 7j+i/2 := 7(2^7+1/2) > 7ij+i/2:= 71 (2^7+1/2) • 

The symbols A± are spatial difference operators: A_Vj := Vj — V,-i and A+l/,- := V^+i — Vj, and we use 
the Engquist-Osher flux [40] 



f{j,u) + f{j,v)- \fu{j,w)\dw 



h{l,v,u) := - 

Note that our pointwise discretization of 7, (3.2), follows the usage of [16, 20, 21, 22], but differs from that 
of [19], where 7 is discretized by cell averages taken over the cells [xj,Xj+i), where Xj := jAx, j £ Z. The 
important point is that in both cases, the discretization of 7 is sto^^ererf with respect to that of the conserved 
quantity u, and this property greatly facilitates the convergence analysis of the numerical schemes. If the 
discretizations were aligned (i.e., not staggered), we would have to deal with more complicated 2x2 Riemann 
problems at cell boundaries. Further discussion of this point is provided e.g. in [ ]. Our particular choice 
of (3.2) (as opposed to forming cell averages) is basically its simplicity. 

The space-time parameters are chosen in such way that we have the following CFL condition (see [20]): 

(3.3) A max |/„(7(a;), u)| + ^ max |A'(u)| < -. 

ue[o,i],a;eR Me[o,i] 2 

which means that At/{AxY must be bounded. On the other hand, when the diffusion term is not considered 
(Example 2 of Section 6), the CFL condition is less restrictive than (3.3), that is 

(3.4) A max \U-f{x),u)\ ^\. 

which means that only At/ Ax must be bounded. 

Let us mention that the scheme also admits a semi-implicit variant, in which the diffusion terms are eval- 
uated at the time level tn+i- This variant has been used for numerical examples in [20], and its convergence 
for a similar equation with a convective flux that does not depend on x, but which is supplemented by 
boundary conditions, has been proved in ['"]. The advantage of a semi- implicit scheme is that it is stable 
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Figure 2. Graded tree data structure. The nodes, leaves and virtual leaves are represented 
by thin, bold, and dotted horizontal lines, respectively. 



under the CFL condition (3.4), which is milder than (3.3), so that much larger time step At could be used. 
However, a semi-implicit version involves the solution of systems of nonlinear equations for each time step, 
and these equations have to be solved iteratively by appropiate linearization. Since we with to keep the basic 
scheme as simple as possible and focus on the multiresolution device, we have decided to avoid this additional 
effort here. Additional complications possibly arise from the fact that we herein implement the scheme on 
an adaptive grid; a semni- implicit variant would, for example, generate nonlinear systems of different size in 
each time step. In general, implicit multiresolution schemes have been explored little so far. 

4. Conservative adaptive multiresolution discretization 

4.1. The graded dynamic tree. The reference standard finite volume scheme described in Section 3 yields 
solutions represented by vectors U" = JJ"'^ containing approximated cell averages on a dyadic uniform grid 
X^ at time i" = nAt. 

An important feature of our scheme is that the differences at different levels, and the solution at different 
levels, are always organized in a tree structure that is dynamic in the following sense: whenever an element 
is included in the tree, all other elements corresponding to the same spatial region in coarser resolutions, are 
also included. The data structure is organized as a dynamic graded tree mainly for sake of data and time 
compression and in particular to be able to navigate through the tree. The adaptive grid corresponds to a 
set of nested dyadic grids generated by refining recursively a given cell depending on the local regularity of 
the solution. The basis of the tree is called root. A node is an element of the tree. In one dimensional space, 
a parent node has two sons, and the sons of the same parent are called brothers. A given node has nearest 
neighbors in each direction, called nearest cousins. A node without sons is called a leaf. For the computation 
of the fiuxes of a leaf, we need s' = 2 nearest cousins in each direction. If these do not exist, we create them 
as virtual leaves. In Figure 2 we illustrate the graded tree structure. 

The nodes of the tree are the control volumes. Following [7], we denote by A the set of indices of existing 
nodes, by C{A) the restriction of A to the leaves, and by A; the restriction of A to a multiresolution level I, 
^ Z s^ L. 

To estimate the cell averages of u on level I from those of the next finer level ^ + 1, we use the projection 
operator Pi^i^i. This operator is exact, unique, and in our one-dimensional case is defined by 



";. 



(Pi 



U, 



i+i^ii-^i+Dj 



1. 



4.2. The multiresolution transform. To estimate the cell averages of a level / -|- 1 from the ones of the 
immediately coarser level I, we use the prediction operator P/^;+i. This operator gives an approximation 
by interpolation of Ui at level ^ + 1. In contrast to the projection operator, there is an infinite number of 
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choice for the definition of P;^;+i, but we impose two constraints: firstly, the prediction is local in the sense 
that the interpolation for a son is made from the cell averages of its parent and the s nearest cousins of its 
parent; and secondly, the prediction is consistent with the projection in the sense that it is conservative with 
respect to the coarse grid cell averages or equivalently, P;+i^/ o P;^;+i = Id. 

For a regular grid structure in one space dimension, we use a polynomial interpolation: 



(4.1a) UI + 1.2J =ULj +^77n{ui,j+m-Uij^m), j ^ I, ■ ■ ■ , N{1) , 

rrt—1 
s 

(4.1b) Ul + i^2j-l = Ul,j ~ 22 lrn{ui,j+in ~ Ulj-m)- 

The order of accuracy of the multiresolution method chosen for our cases is r = 3, which corresponds to 
7i = -1/8 in (4.1a) and (4.1b). 

The detail is the difference between the exact and the predicted value: dij = ui_j — iiij. Given that a 
parent has two sons, only one detail is independent. Then, knowledge of the cell average values of the two 
sons is equivalent to that of the cell average value of the father and the independent detail. Repeating this 
operation recursively on L levels, we get the multiresolution transform on the cell average values M : Ul '~^ 

One of the features of this adaptive multiresolution discretization lies in the possibility to avoid considering 
the prediction error of the numerical flux in the update of the numerical solution, as in Harten's original 
approach. This feature may be seen as an advantage in the frame of equations with discontinuous flux. 

4.3. Multiresolution algorithm. Now we give a brief description of the multiresolution procedure used 
to solve the test problems. 

(1) Initialization of parameters: Model, FV and multiresolution parameters. 

(2) Creating the initial graded tree structure: 

• Create the root of the tree and compute its cell-average value. 

• Split the cell, compute the cell-average values in the sons and compute the corresponding details. 

• Apply the thresholding strategy for the splitting of the new sons. 

• Repeat this until all sons have details below the required tolerance £;. 
DO n — 1 : total _time steps 

(3) Determine the set of leaves and virtual leaves. 

(4) Time evolution with fixed time step: Compute the discretized divergence operator for all the leaves. 
Performing of the space discretization is done in a locally uniform grid (regarding each leaf as a 
control volume of an uniform grid) , so we need only those cell average values which are involved in 
the evaluation of the fluxes for the "edges" of the adaptive mesh formed by the leaves of the tree, 
i.e. we need the leaves and the s' = 2 nearest cousins in each direction. 

(5) Updating the tree structure: 

• Recalculate the values on the nodes and the virtual nodes by projection from the leaves. Com- 
pute the details in the whole tree. If the detail in a node is smaller than the prescribed tolerance, 
then the cell and its brothers are deletable. 

• If some node and all its sons are deletable, and the sons are leaves without virtual sons, then 
delete sons. If this node has no sons and it is not deletable and it is not at level I = L, then 
create sons. 

• Update the values in the new sons by prediction operator from the former leaves. 
END DO n. 

(6) Output: Save mesh, leaves and cell-averages. Deallocate tables and plots. 

With such a process we obtain high order approximation in the smooth regions and mesh refinement near 
discontinuities as a consequence of the polynomial exactness in the multiresolution prediction operator, even 
in the reference finite volume scheme is low order accurate. 
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For a given case of simulation, the performance of the multiresohition method can be assessed by two 
quantities: the data compression rate rj and the speed-up factor V . The data compression rate is defined by 

Nl 

where N^ and |£(A)| are numbers of points of the finest grid and of the leaves in the graded tree, respectively. 
Note that the data compression rate measures the memory compression at a given time of the simulation. 

The speed-up factor is the ratio between the CPU time of the numerical solution obtained by the FV 
method and the CPU time of the numerical solution obtained by the multiresolution method: 

(CPUtime)Fv 

~ (CPUtime)MR' 

5. Error analysis of the adaptive multiresolution scheme 

The main properties of the basic finite volume scheme, i.e., its L^ contractivity, the CFL stability condition 
and the order of approximation in space, allow to derive the optimal choice of the threshold parameter e 
for the adaptive multiresolution scheme. Following the ideas introduced by Cohen et al. [ ] and thereafter 
extended to parabolic equations by Roussel et al. [!)] we decompose the global error between the cell average 
values of the exact solution at the level L, denoted by u^^, and those of the multiresolution computation 
with a maximal level L, denoted by u^r; i'^to two errors 

yO.l) ||Uex "MRII ^ IPcx "FV|| + |PfV "MR||- 

The first error on the right-hand side, called discretization error, is the one of the finite volume scheme on 
the finest grid of level L. It can be bounded by 

(5.2) ||s^^-u^vHC2-"-^, C>0, 

provided that a is the convergence order of the finite volume scheme. The classical approach of Kuznetsov 
[()()] allows to obtain a = 1/2 for a hyperbolic scalar equation. Excepting the discontinuity due to the 
degeneracy, we can anticipate that the value a = 0.5 is a pessimistic estimate of the convergence rate for 
our case. Unfortunately to our knowledge, no theoretical result of convergence rate for numerical schemes 
for strongly degenerated parabolic equations is available so far. Some numerical tests in ["<)] give a « 0.6, 
slightly over our chosen value. 

For the second error, called perturbation error, Cohen et al. [11] assume that the details on a level I are 
deleted when smaller than a prescribed tolerance £;. Under this assumption, they show that if the numerical 
scheme, i.e. the discrete time evolution operator is contractive in the chosen norm, and if the tolerance ei 
at the level / is set to £; — 2'~^e, then the difference between finite volume solution on the fine grid and the 
solution obtained by multiresolution accumulates in time and satisfies 

T 

At" 
where T = nAt and n denotes the number of time steps. 

On the other hand, denoting by |/| the size of the domain and Ax the smallest space step, we have 
Aa; = |/| 2^^. Thus, according to the CFL condition (3.3), the time step At must satisfy 

\J\2 2-2L-1 

At< 



(5.3) W^FV-UunW^C—e, C > 0, 



|/|2-^ max \fuh{x),u)\+ max \A\u)\ 

ue[OA],x£R "£[0,1] 



If we want the perturbation error (5.3) to be of the same order as the discretization error (5.2), we need that 
ex 2^"^. Following Cohen et al. [II], we define the so-called reference tolerance as en := 2^"^Ai. This 



At 

gives 



(5.4) £R = C- 



2-(a + l)l 



I max \fu{j{x),u)\ + 2 ^ max \A\u)\ 

uG[0,l],a;eR Me[0,l] 
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x[mi] 



Figure 3. Example 1 (traffic model): three-dimensional plot of the numerical solution. 



tfinal [h] V 



■q 



L^ error L^ error L°° error 



0.05 6.38 4.5511 5.16x10-^ 6.22x10^5 5.64x10"^ 

0.10 6.99 4.2140 4.57x10"'' 2.41x10"^ 8.16x10"'^ 

0.15 7.84 7.8168 7.21x10"'' 5.12x10"^ 7.23x10"'* 

0.20 9.01 17.3559 1.14x10"^ 2.47x10"^* 3.86x10"^^ 

Table 1. Example 1 (traffic model): Corresponding simulated time, speed-up factor V , 
compression rate, and normalized errors. L = 10 multiresolution levels. 



For the case A{u) — (see Example 2 of Section 6), the reference tolerance must be taken as 

^ -1 



max \fuh{x),u) 



because of the less restrictive CFL condition. To choose an acceptable value for the factor C, a series of 
computations with different tolerances are necessary. 



6. Numerical results 

6.1. Example 1: Diffusively corrected kinematic traffic model with changing road surface con- 
dition. Our numerical example for this model has been chosen in such a way that results can be compared 
with simulations shown in Example 5 of [ ]. The velocity function is given by (2.2) with Umax = 220 cars/mi, 
C — e/7 = 0.38833 and Wmax = 70mph, so that 



Vina.xU for < u < Uc = exp(— 1/C)m„ 

(6.1) f{u) = \ Umax(e/7)uln(Uniax/u) for Uc <U< Umax, 

otherwise. 



16.7512 cars /mi. 
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(a) 



(b) 




.5 2 2.5 3 



-0.5 0.5 



1.5 2 2.5 3 



(c) 



(d) 





.5 2 2.5 3 



Figure 4. Example 1 (traffic model): (a, c) numerical solution and (b, d) positions of the 
leaves at (a, b) i = 0.05 h and (c, d) t = 0.1 h. 



We choose fH,. 



70 mph and v^^^ 
Po{x) 



25 mph. The initial density is chosen as 

J 100 cars/mi for —2 mi < x < — Imi, 
I otherwise. 



The integrated diffusion coefficient A{u) resulting from our choice of parameters satisfies A{u) = for 
< u < Mc = 16.7512 cars/mi, and has an explicit algebraic representation [IG, 58]. 

In Example 1, we consider an initial convoy of cars traveling on an empty road, and wish to see how the 
convoy passes through the reduced speed road segment. The numerical solution obtained by our method is 
represented in a three-dimensional plot in Figure 3 and shown at four different times in Figures 4 and 5. 
These figures also display the corresponding position of the leaves. For these four times, Table 1 displays the 
corresponding values of the speed-up factor V, the compression rate t], and normalized approximate errors. 
These errors and the speed-up factor are measured with respect to a fine grid calculation (no multiresolution) 
with Nl — 2^^ cells. (We further comment on the behaviour of V and 77 in the discussion of Example 3.) 

For this example, we take an initial dynamic graded tree, allowing L = 10 multiresolution levels. We use 
a fixed time step determined by A = 0.0003 h/mi, thus At = Xh^. The prescribed tolerance er is obtained 
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Figure 5. Example 1 (traffic model): (a, c) numerical solution and (b, d) positions of the 
leaves at (a, b) i = 0.15h and (c, d) t = 0.2 h. 



from (5.4), where the constant C for this example corresponds to a factor C — 10^, so e = 0.301 and the 
thresholding strategy is Ek — 2'^^^e. 

The errors in L^ norm between the numerical solution obtained by our multiresolution scheme for different 
multiresolution levels L, and the numerical solution by finite volume approximation in a uniform fine grid with 
2^^ control volumes, are depicted in Fig. 6. In practice, we compute the error between the numerical solution 
obtained by multiresolution and the projection of the the numerical solution by finite volume approximation. 
We also observe the same slope (= 0.8819) between finite volume and multiresolution computation in the 
L^ error of Figure 6. 

6.2. Example 2: Clarifier-thickener treating an ideal suspension {A = 0). For Example 2, we 
choose the same parameters as in [)(), -If], so that results can be compared. In particular, we consider an 
ideal suspension that does not form compressible sediments, i.e., we set A = 0, so that the model considered 
in this example actually corresponds to the first-order equation (1.3). 

We consider a clarifier-thickener unit that is initially full of water by setting uo(x) =0. At i = 0, we 
start to fill up the device with feed suspension of concentration u-p = 0.8. We also consider .tl = —1 and 
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Figure 6. Example 1 (traffic model): L^ errors. 




Figure 7. Example 2 (clarifier-thickener model with ^ = 0): three-dimensional plot of the 
numerical solution. 



a;R = 1 and we assume that the mixture leaving the unit at xl and xr, is transported away at the bulk flow 
velocities (Zl = — 1 and g^ = 0.6. The suspension is characterized by the function J{u) given by (2.9) with 

Woo = 27/4, C = 2 and Umax = 1- 

We use an initially graded tree with i = 10 multiresolution levels and a reference tolerance of e = 
4.15 X 10^'^. The finest grid has Nl — 512 control volumes and we choose a factor A = 1/16. Observe 
that the visual grid used to display Figure 7 coincides with the computational grid in x direction, but in t 
direction, only every 50th profile is plotted. 
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Figure 8. Example 2 (clarifier-thickener model with A = 0): numerical solution (a, c) and 
position of the leaves (b, d) at t = 1 and t = 2. 



ifinal V 



V 



L^ error L^ error L°° error 



1 8.42 6.4362 2.47x10-* 6.31x10-'' 8.49x10"^ 

2 9.36 8.0315 4.11x10-* 8.47x10"'' 2.40x10-" 

3 10.21 8.7850 3.42x10-* 1.84x10^3 6.74x10-* 

4 10.94 8.7850 4.18x10-* 1.10x10"^ 1.26x10'^ 

Table 2. Example 2 (clarifier-thickener model with A = 0): Corresponding simulated time, 
speed-up factor V, compression rate, and normalized errors. L = 10 multiresolution levels. 



For Example 2, we use as a reference solution a fine grid computation with 2^^^ control volumes. Table 2 
lists the behaviour of the error and the gain in computational effort and data storage for different times. 
Also, analogously to Example 1, we can observe in Figure 10 that the plots of the L^ error, which is measured 
here for i = 2, have the same slopes. 



16 



BURGER, RUIZ, SCHNEIDER, AND SEPULVEDA 



(a) 



(b) 



,tm*^..00^^^^^mmmmmm*mm^ 



-0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 



-0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 



(c) 



(d) 













,iiffl***»*«**«****^*********'^ 


0.8 




* * *^Hf*rf 


0.6 




- 


0.4 




- 


0.2 




^ 




":'':''':'':'':''':'':'':''':'':'':i 



-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 




-0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 



Figure 9. Example 2 (clarifier-thickener model with A = 0): (a, c) numerical solution and 
(b, d) positions of the leaves at (a, h) t = 3 and (c, d) t = 4. 



ifinal [s] 


V 


V 


L^ error 


L^ error 


L°° error 


10000 


7.88 


4.1787 


3.67x10-4 


8.41x10-5 


6.73x10-4 


25000 


9.01 


4.4265 


4.82x10-4 


9.32x10-5 


8.29x10-4 


50000 


10.74 


4.4734 


6.30x10-4 


1.24x10-4 


1.07x10-3 



Table 3. Example 3 (clarifier-thickener model with A ^ 0): Corresponding simulated time, 
speed-up factor V , compression rate, and normalized errors. L = 9 multiresolution levels. 



6.3. Example 3: Clarifier-thickener treating a fiocculated suspension {A ^ 0). The parameter of 
the flux is the same as in Example 2, and the function adu) is given by (2.12) with ctq — l.OPa, Uc = 0.1 and 
/3 = 6. The remaining parameters are Ag = 1660 kg/m"^ and g — 9.81 m/s^. Note that for (2.9) with /3 G N, 
the function A(u) has an explicit closed-form representation, see [(>]]. The reference numerical scheme is 
(3.1) with A = 40s/m. 
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Figure 10. Example 2 (clarifier-thickener model with A = 0): L^ errors. 
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Figure 11. Example 3 (clarifier-thickener treating a flocculated suspension): two views of 
the time-space representation of the numerical solution. 



Our simulation corresponds to the choice q^ — 2.5 x 10^^ m/s and ql = —1.0 x 10^'''m/s. The feed 
concentration corresponds to up — 0.086. Figure 11 shows the numerical solution until t = 50000 s. In this 
case we consider the device with an initial concentration distribution of uq(x) — u^, x ^ [a^Lj^;^] and we 
can observe the initial stage of the fill-up process. For this example, we take an initial dynamic graded tree, 
allowing L — 8 multiresolution levels, and for the reference tolerance we use C = 10"'^, so £r = 2.24 x 10""*. 

For Example 3 we use as a reference solution a fine grid computation with 2^^ control volumes. Table 3 
again displays the behaviour of the error and the gain in computational effort and data storage for different 
times. Also, analogously to the previous examples, we observe in Figure 14 the same slope between the L^ 
errors for the finite volume and multiresolution methods. This error is masured here at i = 25000 s. 

Note that in all numerical examples, the speed-up factor V increases as ifinai is increased, even if the 
data compression rate rj remains constant, which approximately is the case in Table 3, or even decreases. 
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Figure 12. Example 3 (clarifier-thickener treating a flocculated suspension): (a, c) numer- 
ical solution (stars) and initial condition (dashed), (b, d) positions of the leaves (plus) at (a, 
h) t = 10000 s, (c, d) i = 25000 s for the transition to a steady state from uq = 0.09. The 
horizontal solid line in (a, c) denotes the critical concentration Uc = 0.1. 

as we see, for example, by comparing the values of r] for ifinai — 0.05 h and ifinai = 0.10 h in Table 1. The 
explanation of this discrepancy is that while ij measures the quality of performance of the multiresolution 
seen at the instant t = ifinai, the speed-up factor V is referred to the total time of simulation and also 
includes the "overhead" required by initializing the graded tree in step (2) of the multiresolution algorithm. 
The initialization requires a flxed amount of CPU time, which is independent of the number of total time 
steps (which is porportional to tfinai, since we consider At to be fixed). On the other hand, a standard FV 
method on a fixed grid will always require CPU time proportional to the number of time steps. This explains 
why even if rj does not change significantly, we observe an inprovement of the speed-up factor V as ignai is 
increased. 



7. Conclusions 

Before discussing our results, we comment that the standings of both applicative models are slightly 
different. Numerous mathematical models have been proposed for onc-dircctional flows of vehicular traffic; 
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Figure 13. Example 3 (clarifier-thickener treating a flocculated suspension): (a) numerical 
solution (stars) and initial condition (dashed) , (b) positions of the leaves (plus) at t = 50000 s 
for the transition to a steady state from uq = 0.09. The horizontal solid line in (a) denotes 
the critical concentration Uc = 0.1. 
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Figure 14. Example 3 (clarifier-thickener model with ^ ^ 0): L^ errors. 



reviews of this topic are given in the monographs by Helbing [(i2], Kerner [(i.H] and Garavello and Piccoli 
[G4], as well as in the articles by Bellomo et al. [()5, 66, 67]. These and other works vividly illustrate that 
the number of balance equations (for the car density, velocity, and possibly other flow variables) that form 
a time-dependent model based on partial differential equations, as well as the algebraic structure of these 
equations, is a topic of current research. Fortunately, all these models are spatially one-dimensional, and a 
circular road with periodic boundary conditions provides a setup that is both physically meaningful (since 
the flow is horizontal) and easy to implement for numerical simulation. This setup, on one hand, is widely 
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used to compare different traffic models, and on the other hand allows to assess the local influence and 
long-term behaviour of nonlinearities and inhomogeneities such as the ones introduced in Section 6.1. 

While the traffic problem highlights the use of the scheme used herein to explore different models, the 
clarifier-thickener model calls for an efficient tool to perform simulations, on one hand, related to clarifier- 
thickener design and control [20, G(S], and on the other hand, to parameter identification calculations [G9, 70]. 
In fact, depending on the parameters, clarifier-thickener operations such as fill-up may extend over weeks 
and months [68], and require large simulation times, while the parameter identification procedures in [69, 70] 
proceed by solution of an adjoint problem, which needs storage of the complete solution of the previously 
solved direct problem. Clearly, methods that imply savings in both computational time and memory storage, 
such as the multiresolution scheme presented herein, are of significant practical interest for the clarifier- 
thickener model. 

Both mathematical models considered herein exhibit three types of fronts that typically occur in solutions 
of (1.2), namely standard shocks (i.e., discontinuities between solution values for both of which (1.2) is 
hyperbolic), hyperbolic-parabolic type-change interfaces (such as the sediment level in Example 3), and 
stationary discontinuities located at the discontinuities of 'y{x). The basic motivation for applying a finite 
volume multiresolution scheme is that this device is is sufficiently flexible to produce the refinement necessary 
to properly capture all these discontinuities, and leads to considerable gains in storage as can be seen from the 
sparsity of the graded trees in our numerical examples. Moreover, Figure 6 confirms that we may effectively 
control the perturbation error, in the sense that the error of the resulting finite volume multiresolution 
scheme remains of the same order as that of the finite volume scheme on a uniform grid. We recall from 
Section 5 that the feasibility of this control depends on an estimate of the convergence rate of the basic 
discretization on a uniform grid, which is an open problem for strongly degenerate parabolic equations. 

Although our numerical results look promising, they still alert to some shortcomings that call for im- 
provement. The most obvious one is the limitation of the time step according to the spatial step size of 
the finest grid, which can possibly be removed by using a space-time adaptive scheme such as the recent 
finite volume multiresolution schemes by Stiriba and Miiller [ ]. On the other hand, the basic finite volume 
scheme accurately resolves the discontinuities of the solution sitting at the jumps of 7(3;) at any level of 
discretization; these discontinuities are not approximated by smeared transitions (as are discontinuities at 
positions where 7(x) is smooth), see [1!)]. This means that the refinement the multiresolution produces 
near these discontinuities, which is visible in Figures 8 and 9, and which is based on the adaptation of the 
refinement according to features of the solution (but not of 7(2;)), is possibly unnecessary, and that a more 
efficient version of the present method may be feasible. 
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